Computing the merger of black-hole binaries: the IBBH problem 
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Gravitational radiation arising from the inspiral and merger of binary black holes (BBH's) is a 
promising candidate for detection by kilometer-scale interferometric gravitational wave observato- 
ries. This paper discusses a serious obstacle to searches for such radiation and to the interpretation 
of any observed waves: the inability of current computational techniques to evolve a BBH through 
its last ~ 10 orbits of inspiral (~ 100 radians of gravitational- wave phase). A new set of numerical- 
relativity techniques is proposed for solving this "Intermediate Binary Black Hole" (IBBH) problem: 
(i) numerical evolutions performed in coordinates co-rotating with the BBH, in which the metric 
coefficients evolve on the long timescale of inspiral, and (ii) techniques for mathematically freezing 
out gravitational degrees of freedom that are not excited by the waves. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.-s 
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I. MOTIVATION 

Among all gravitational wave sources that theorists 
have considered, the one most likely to be detected first 
is the final inspiral and merger of binary black holes 
(BBH's) with masses M x ~ M 2 ~ 1O-2OM @. Detailed 
analyses of the evolution of stellar and black-hole popu- 
lations predict event rates as high as ~ one per year 
in the first LIGO/ VIRGO interferometers (2002-2003) 
and a thousand per year in enhanced interferometers for 
which research and development is currently under way, 
but the rates could also be far lower than this. 

Optimal search techniques require prior information 
about the gravitational waveforms. The waveforms from 
the early binary inspiral phase, when the holes are far 
apart, are calculated by a post-Newtonian (PN) expan- 
sion. The merger phase, beginning at the innermost sta- 
ble circular orbit, will be calculated by numerical relativ- 
ity. Unfortunately, there is a gap || between the failure 
of the PN expansion (which, for concreteness, we take 
to occur when its Taylor series makes a 2% error in the 
energy loss rate Q) and the beginning of merger. Fill- 
ing this gap is called the Intermediate Binary Black Hole 
(IBBH) Problem Q. 

We estimate ff the PN failure point, for calcula- 
tions at 3PN order [0(w 6 ) beyond Newtonian gravity 
and quadrupolar radiation reaction] , to be at the orbital 
speed v = (Affi) 1 / 3 ~ 0.3 (where M is the system's total 
mass, is its orbital angular velocity, and G = c = 1); 
there the remaining time to merger, remaining number 
of orbits, and remaining number of gravitational-wave 
radians are T ~ 1200M , N OThits ^ 8, and $ ~ 100. For 
2.5PN calculations, the PN failure is at v ~ 0.25 where 
T ~ 5000M, Ambits ^ 20, and $ ~ 250. For optimal 
detection of the waves, the waveform must be accurately 
modeled in the IBBH gap @. The wave frequency in 
this gap is / = tt/ir ~ (50 to 200 Hz)(2OM /M), which 
is the band of optimal LIGO/VIRGO sensitivity. This 



adds urgency to the IBBH problem. 

For numerical simulations of the merger phase, the 
conventional approach uses asymptotically inertial coor- 
dinates in which the dynamical timescale, Td yn ~ M , is 
set by the task of moving the holes across the coordinate 
grid. It is unlikely that, in the next several years, this ap- 
proach will be able to evolve a BBH through the gap for 
the required > 1200 dynamical time scales. This moti- 
vates exploring alternative procedures for computing the 
evolution and waves during the IBBH phase. 

One possible method of extending the PN expansion 
into and through the IBBH region is to augment it with 
Pade approximants || . Preliminary results || , in which 
PN and exact waveforms in the test-mass limit are com- 
pared, give cause for optimism that the waveforms' phase 
evolution can be computed with adequate accuracy via 
PN Pade-approximants all the way through the IBBH 
region. However, there is little hope, via PN Pade ap- 
proximants, to evolve the binary's internal spacetime ge- 
ometry in the IBBH region and thereby provide (i) initial 
data for numerical relativity's analysis of the merger, and 
(ii) a connection between those initial data and the bi- 
nary's early inspiral properties (masses, spins, orbit). For 
these crucial issues we must turn elsewhere. 

In this paper we explore an alternative strategy [gj : nu- 
merical relativity computations performed not in asymp- 
totically inertial coordinates (as is normally done), but 
instead using spatial coordinates which co-rotate with 
the holes' orbital motion and a temporal slicing which 
adjusts, as the potential well between the holes deepens, 
so as to keep all the metric coefficients as slowly evolv- 
ing as possible. In such coordinates it is reasonable to 
hope to achieve a timescale Td yn for dynamical evolution 
of the metric coefficients that is of order the timescale t* 
on which radiation reaction drives the holes together.^] 



"If the holes are spinning with axes inclined to the orbital 
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Since the orbital frequency changes by only a factor of 
~ 2-3 through the IBBH phase, this phase may last only 
~ 3 dynamical timescales in the co-rotating frame — an 
enormous reduction from the > 1200 timescales in the 
asymptotically inertial coordintcs of standard numerical 
relativity. 

Although the metric coefficients' true dynamical 
timescale will be t* in these co-rotating coordinates, 
numerical approximations may excite spurious gravita- 
tional waves with wavelengths of order the spatial grid 
size. Some numerical schemes will be forced to take time 
steps shorter than the grid size in order to suppress these 
modes and to keep the numerical evolution stable; cf. 
the Courant condition. A good numerical scheme should 
freeze out these unphysical modes and stabilize the evo- 
lution while using long time steps. Correspondingly, a 
concrete implementation of our strategy must include 
two elements: first, a method to choose the lapse and 
shift so the coordinates co-rotate with the binary; sec- 
ond, a numerical scheme that evolves stably with time 
steps constrained only by . Such a scheme differs from 
that of previous co-rotating neutron-star-binary calcula- 
tions ||, which have not evolved the gravitational field 
but instead computed sequences of equilibria. 



dt- This is a vague statement, which we make precise 
by thinking of the left side of Eq. (||) as a velocity, con- 
structing a kinetic energy from this velocity, and choosing 
a lapse and shift that minimize this kinetic energy We 
will discuss several such action principles for a and f3j in 
the next two subsections. 



A. The Minimal-Strain Lapse and Shift 

In the spirit of Smarr and York's minimal distortion 
shift, we can construct an action principle based on min- 
imizing the Lie derivative of the spatial metric jij rather 
than the spacetime metric. Specifically, we presume that 
the numerical evolution has proceeded up to some slice of 
constant time t that has intrinsic metric 7^ and extrinsic 
curvature , and we choose the lapse a and shift 0i on 
this slice so as to minimize the positive definite action 



h[a,/3k] = J d 3 V7 



ik jl • 

Jijl T Ikl 



(3) 



Here jij — djij/dt (the Lie derivative of jij along dt) 
must be expressed in terms of Kij , a, f3j via the standard 
relation 



II. CHOOSING THE LAPSE AND SHIFT 

Numerical relativity is based on a 3 + 1 decomposition 
of the metric: 



ds 2 = 



t 2 dt 2 + 7y (dx i + f3 l dt)(dx° + (3 3 dt) . (1) 



Here a is the lapse function, j3 l is the shift vector, and 
7y is the intrinsic metric of the 3-dimensional slices of 
constant time t. The lapse and shift are specified freely 
during the evolution, thereby fixing the spacetime coor- 
dinates. 

We propose to construct the initial IBBH co-rotating 
coordinates and metric from the PN metric near the PN 
failure point by adjusting the lapse a and shift f3j so 
as to make the metric coefficients evolve on the inspiral 
timescale r*. Subsequently a and j3j must be chosen so 
as to make the coordinate time derivatives of all the met- 
ric coefficients stay small — or, equivalently in coordinate- 
independent language, to make 



C 9t g =s , 



(2) 



where Cg t g is the Lie derivative of the spacetime metric 
g with respect to the coordinate system's time generator 



angular momentum, then in these coordinates the evolution 
timescale may be shorter: r* ~ (spins' precession period). For 
simplicity we shall ignore this possibility, though our analysis 
presumably can be adapted to handle it. 



jij = -2aKij + 2D(if3j) , 



(4) 



where Di is the spatial gradient compatible with the 3- 
metric 7y. By minimizing the resulting action with re- 
spect to variations of a and /3,-, we obtain four coupled 
equations: 



K ij [-2aKij +2Di(3j] = 
Di[-2aK i j+2D {i f3j ) ] = 



(5a) 
(5b) 



Equation (^aj) is easily solved to give a in terms of (3j. 
When that a is inserted into Eq. (|5b|), the result is a 



linear, homogeneous differential equation for f3j. If the 



lapse were not fixed via Eq. (5a) but instead were chosen 
independent of /3j, e.g., via maximal slicing, then the shift 
equation (5b) would reduce to the minimal strain shift 
of Smarr and York 0. We therefore refer to Eqs. (||) as 
minimal strain equations. 

Notice the geometrical nature of the spatial coordi- 
nates carried by this lapse and shift: The action prin- 
ciple ([|) minimizes the rate of change, along dt, of the 
infinitesimal proper distance between neighboring points 
at fixed spatial coordinates. This (presumably) will be 
achieved, in the binary itself, by making the coordinates 
co-rotate with the holes, and in the radiation zone by at- 
taching the spatial coordinates to the wave pattern, i.e., 
by (almost) freezing the wave pattern into the spatial 
coordinate grid. A direct consequence is that evolution 
along dt is nearly shape and volume preserving. 

In the IBBH problem, this approach is not without 
shortcomings: there is no guarantee that the minimal 
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strain equations, which are solved on each spatial slice, 
will force the lapse and shift to evolve slowly. On the 
other hand, if the initial data are constructed in coordi- 
nates that are close to co-rotating (as they will be using 
the known PN metric), and if appropriate slow-change 
boundary conditions are enforced on /3j near the holes' 
apparent horizons and at the outer edge of the coordinate 
grid, then it is reasonable to expect a and to evolve on 
the same slow timescale t» as the spatial metric %j . This 
is because a and (3i inherit their dynamics from the time 
evolution of 7^ and 2£y. Note that the minimal strain 
equations become degenerate for time-symmetric initial 
data; such a situation will not arise in the IBBH prob- 
lem. A method of enforcing a variant of Eq. (|5a| ) where 
K % i is replaced by 7^ has been explored by Balakrishna 
et al. § 

The following (far from rigorous) argument makes it 
seem likely that this scheme will succeed for the IBBH 
problem. The IBBH spacetime has an "almost Killing 
vector field" £, which embodies co-rotation and satisfies 



C € g = s ~ A/7 



(0) 



Here A ~ M is the length scale over which the spacetime 
curvature varies, and r* ^> A is the inspiral timescale, 
so s is small. In terms of the 3+1 spacetime foliation 
being generated by the minimal-strain lapse and shift, 
we can decompose £ into a spatial piece B and a piece 
in the direction h normal to the surfaces of constant t: 
£ = An + B, where B ■ n — by definition. We wish to 
determine the effectiveness of the minimal strain equa- 
tions at attaching the coordinate grid to £, i.e., at mak- 
ing £ = dt or equivalently A = a and B = (3. First, we 
project Eq. (0) into the spatial slice S to get 



ID^Bj) - 2AK l3 = Sij . 
Next, because a and (3 l satisfy Eqs. (||), 
TP [—2KijK kl Dkf3i / K mn K mn + 2D (i i) ] 



(7) 



(8) 



Finally, substituting (3 l = B % — b % into this equation and 
using Eq. (Q) we find that b l , the difference between the 
minimal-strain shift and the shift we would like, satisfies 



D 3 [-2K l3 K kl D k bi/K mn K mn + 2D (l b 3) ] 
= D^[-K ij (K kl s kl /K mn K mn ] 



(9) 



Assuming (without proof) that the boundary value 
problem for Eq. (||) is well posed, we see that there ex- 
ists a solution to Eq. (||) that will scale as b 1 ~ A/t„; 
Eq. ( |5a| ) then implies that a — A ~ A/t*. Therefore, the 
minimal-strain shift and lapse can make dt equal to the 
almost Killing vector field £ that embodies co-rotation, 
aside from fractional differences of order A/r*, as we de- 
sired. 

Notice that if £ = Ah + B is a Killing vector field on 
the spacetime then = 0, and b % — is a trivial solution 
to Eqs. (0) corresponding to a = A and (3 l — B 1 . 



B. Other Choices of Lapse and Shift 

There is much freedom in choosing the lapse and shift 
to achieve the goal of slowly evolving metric coefficients. 
Another class of action principles that might work is 
based on minimizing an integral over spacetime rather 
than over 3-space as in Eq. (||) . Let dt — an + (3 be the 
vector field to which our coordinates are tied, and de- 
note the Lie derivative of the 4-metric along dt by — 
CQ t g fjLl> . Let v be some other vector field independent of 
dt, and from it construct the tensor — + v^v u . 
Then our class of actions is 



h[df,v]= I U^H^H^j pa ) 

J M 



(10) 



On varying dt , while holding v and the spacetime metric 
fixed, we arrive at 



(11) 



This is a dynamical system of equations for the lapse and 
shift. Certain values of v might be considered most nat- 
ural. If v = \/2x (some unit timelike vector), then H£ u is 
positive definite and there is a solution of Eqs. (0) that 
truly minimizes the action. If v = 0, then Eqs. (|ll|) are 
a simple conservation law, but the action is not positive 
definite. It is trivial to show that a spacetime Killing 
vector field is a solution to Eqs. ( pd| ) independent of the 
choice of v, and straightforward to extend the analysis 
of Eqs. (0)~@ to show that for the IBBH problem one 



of the solutions of ( 
vector field" of Eq. 



(nl) differs from the "almost Killing 
J) by an amount that scales as A/r*. 
However, neither here nor for our minimal-strain equa- 
tions have we managed to demonstrate that the desired 
solution for dt is an attractor; this needs further study. 



III. NUMERICAL EVOLUTION 

To fully solve the IBBH problem will require combining 
one of our lapse/shift differential equations with the Ein- 
stein equations in some concrete numerical scheme. As 
noted in Sec. |, although the binary's metric coefficients 
should evolve on the long timescale t* in our proposed 
co-rotating coordinate system, there is danger that the 
time steps will be driven down to less than the size of the 
spatial grid by the numerical scheme's attempt to follow 
spurious gravitational waves and/or to control numeri- 
cal instabilities (the Courant condition). To avoid these 
pitfalls while taking time steps controlled only by the in- 
spiral timescale r*, it will be necessary to stabilize the 
integration scheme and freeze out the degrees of freedom 
that are physically present but unphysically excited. 
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A. Freezing Out Unwanted Degrees of Freedom 

It is well known that implicit differencing schemes 
freeze small-scale structure and produce unconditionally 
stable evolution. For this reason, we envisage using an 
implicit scheme to evolve the Einstein equations in co- 
rotating coordinates. 

While implicit differencing may be sufficient to achieve 
a stable evolution on the timescale r*, we propose an 
additional technique that should also help. The idea is 
to convert the ADM equations for 7^- and Kij into a 
parabolic system, thereby removing all spurious waves 
while keeping the real ones (which are nearly frozen into 
the co- rotating coordinates). 

The evolution equations written in the usual ADM 
form are 



-2aK, 



(12a) 



= —DiDja + a[ s Rij + j kl (K l3 K H - 2K lk K l3 )] 



+CpK, 



i] 1 



(12b) 



where 3 i?y is the Ricci tensor constructed from 7^ . This 
first-order system can be r e-ex pressed as a second or- 
der system by solving Eq. ( 12a ) f or th e extrinsic curva- 
ture and substituting it into Eq. 



12b). Since the fields 



evolve on the very long timescale r* in the co-rotating 



frame, and since (dt) 



incrtial 



(dt) 



co— rotating 



VLd,h, the 



terms with two time derivatives in co-rotating coordi- 
nates will be smaller, by a factor ~ l/(f2r»), than at least 
some of those with a single time derivative. Thus, the 
double-time-derivative terms can be neglected (or back- 
differenced, if desired, so they are treated as sources aris- 
ing from data on previous time slices). In particular, the 
term involving 7^ can be neglected (or back-differenced). 
(Since a and /3j are not dynamical fields, we suggest 
that their time derivatives also be back-differenced.) The 
resulting parabolic system of equations for 7^ can be 
evolved using an implicit scheme, which should be stable 
for large time steps. 



frame. Such matching could conceivably be done around 
each of the holes and at an outer boundary in the radi- 
ation zone jjuj. It may be possible to impose outgoing- 
wave boundary conditions as a constraint on the spatial 
derivatives of 7^ at the outer boundary. The shift there 
is (3 ~ ^ld r f >1 where 8$ is the generator of rotations in 
the orbital plane, and outgoing waves are constant along 
d r + (3 (aside from their 1/r amplitude falloff), whereas 
spurious ingoing waves would be constant along d r — {3. 

We also need boundary conditions for the differential 
equation used to compute the shift. Fortunately, these 
seem easy to construct. Since the matching to charac- 
teristics would be done on the history of closed spatial 
2-surfaces on the outer boundary and around each black 
hole, the derivations of the equations for the lapse and 
shift can be repeated for these surfaces. For example, 
the minimal strain equations, obtained by the variations 
of an action of the squared velocity of the metric on the 
2-surface [the analog of Eq. (||)] , are 



k ab [-2ak ab + 2V a (3 b ] = 
V a [-2ak ab + 2V {a /3 b) } =0 



(13a) 
(13b) 



where k ab is the extrinsic curvature of the 2-surface em- 
bedded in its history, T> a is the covariant derivative op- 
erator compatible with the metric on the 2-surface, and 
P a is the 3-dimensional shift vector projected into the 
2-surfacc. Since the boundary is closed, these equations 
can be solved to obtain the lapse function and the tan- 
gential components of the shift vector on the boundary. 
The component of the shift that is normal to the bound- 
ary would be set to zero, so the evolution will not at- 
tempt to follow wave crests as they leave the numerical 
grid. Among the various solutions to Eqs. Jl5|), the one 
we want will be that which is closest to the solution on 
the previous time slice. This boundary solution, com- 
bined with our 3-dimcnsional differential equations for 
the lapse and shift, presumably will produce the desired 
a and which evolve on the slow time scale r* and 
co-rotate with the holes. 



B. Initial Data and Boundary Conditions 

To solve the IBBH problem, we must specify suitable 
initial data and boundary conditions in addition to for- 
mulating an evolution scheme and a method of choos- 
ing the lapse and shift. One can construct the initial 
data, just before the PN failure point, by using matched 
asymptotic expansions to join the post-Newtonian exte- 
rior metric onto the metrics of two tidally distorted Kerr 
black holes, and by then transforming to co-rotating co- 
ordinates ||. 

The method of Cauchy characteristic matching (e.g., 
Ref. JTof) seems a promising candidate for construct- 
ing boundary data for evolution of 7^ in the co-rotating 



IV. WORRIES 

Two conceivably serious difficulties with our approach 
are: (i) in our co-rotating reference frame, the almost 
Killing vector becomes spacelike beyond the speed-of- 
light surface, which might cause problems for the evolu- 
tion; (ii) when the second time derivatives are discarded, 
the resulting evolution might not represent the true evo- 
lution of the spacetime. We doubt, however, that these 
difficulties will actually arise, since we have seen no sign 
of them in a toy problem that retains the relevant fea- 
tures of the IBBH problem. 

Our toy problem is a rotating, radiating sphere of 
scalar charge in flat spacetime. We begin in standard 
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spherical polar coordinates so the line element is given 
by Eq. (Q) with x % = {r,6,<p}, 7^ = diag[l, r 2 , r 2 sin 2 9] 
and (3 l = 0. The sphere has radius R, angular veloc- 
ity f2(f), moment of inertia I and scalar charge density 
p[r, 9,(t>-<f>(t)} where $(t) = /' dt'f2(f ). The scalar field 
\1> produced by this charge distribution satisfies the wave 
equation 

g al3 W a W^ = 47T P e(R-r) , (14) 

where Q(R — r) is the step function which is unity in- 
side the sphere and zero outside. As the sphere rotates, 
scalar waves arc radiated to infinity, decreasing its an- 
gular momentum IQ according to the radiation-reaction 
equation 

1^=4*/ d 3 x^f [ P (d^) Q(R - r)] . (15) 

For a quadrupolar distribution of the scalar charge, this 
model can be reduced, via separation of variables, to a 
(l+l)-dimensional problem, which we have evolved in 
the inertial frame by standard finite-difference methods 
for hyperbolic systems and time steps constrained by the 
Courant condition. 

We have also evolved this 1 + 1 system in the co- 
rotating frame {r,9,cf> — <f> — $(£),£} with $(f) inferred 
from the radiation-reaction equation ( |i~5| ) and not from 
some variant of our minimal-strain equations. In this 
co-rotating evolution, we discarded the (small) second 
time derivatives of \& and applied an implicit differenc- 
ing scheme to the resulting parabolic system. We suc- 
ceeded to make the time steps as large as the rotational 
timescale of the charged sphere (much larger than the 
Courant condition would allow), and we believe that the 
time steps could be made as large as the radiation reac- 
tion timescale (which is much larger than the rotational 
timescale) , since the limiting factor was the simple outer 
boundary condition that we used. There were no numer- 
ical instabilities. Moreover, no numerical problems were 
encountered at the speed-of-light surface r = (fisinfl) -1 
[and we expected none, since the transformation to co- 
rotating coordinates does not change the fact that the 
wave equation (|l4| ) is manifestly hyperbolic]. There was 
good agreement between the results computed in the in- 
ertial and co-rotating frames. Further details will appear 
elsewhere G|. 

Additional evidence that evolution in rotating coor- 
dinates need not cause problems comes from the work 
of Bishop et al. jo. They have used Schwarzschild 
spacetime transformed into rigidly rotating coordinates 
as a test bed for their numerical evolutions of the Ein- 
stein equations, and they encountered no problems at the 
speed-of-light surface. They report that their character- 
istic code "can routinely handle such superluminal gauge 
flows" at large distances from the Schwarzschild black 
hole, even though the coordinate grid moves slower than 
the speed of light near the hole. 



Based on these results, it seems likely that an imple- 
mentation of the methods presented here will allow a nu- 
merical evolution of a binary system through the IBBH 
phase. Since a better understanding of this phase is 
important — and perhaps critical — for the LIGO /VIRGO 
detection of waves from binary black hole systems, 
and since such systems are highly promising candidate 
sources for LIGO and VIRGO, we hope to inspire re- 
searchers in numerical relativity to address the IBBH 
problem via our proposed techniques or others. 
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